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NEAR-NEIGHBOR  ALGORITHMS  FOR  PROCESSLNG  BEARING  DATA 


I.  Introduction 

One  of  the  more  pressing  problems  in  Correlator-Tracker  (CT)  algorithms  is  the  ef¬ 
ficient  use  of  sensor  data  consisting  of  only  the  bearings  of  objects  relative  to  the  sensor. 
By  assumption,  no  range  information  is  available,  or  the  range  is  much  less  accurately 
known  than  the  bearing.  Infrared  sensors  will  produce  bearing  data,  for  example.  Here 
we  present  recent  progress  in  processing  and  using  such  data.  The  paper  emphasizes  use 
of  data  structures  designed  for  ready  access  of  information  on  the  near  neighbors  of  each 
object,  such  as  a  track  or  contact  report,  within  an  ensemble.  Implicit  in  the  discussion  is 
the  perceived  need  for  optimization  on  highly  parallel  or  vector  computers  of  the  present 
and  future  in  order  to  handle  tens  or  hundreds  of  thousands  of  objects  as  rapidly  as  possi¬ 
ble.  The  goal  is  to  perform  the  same  calculation  simultaneously  on  large  subsets  of  these 
objects. 

Since  the  likelihood  of  association  of  a  track  with  a  sensor  observation  or  “report-’  de¬ 
creases  with  track-report  separation,  one  usually  assumes  that  rapid  access  of  near  neigh¬ 
bors  will  drastically  reduce  the  work  required  in  updating  tracks  through  the  use  of  new 
sensor  data.  The  next  section  discusses  the  conditions  under  which  the  use  of  near-neighbor 
(NN)  algorithms  is  actually  advantageous.  This  discussion  also  indicates  conditions  under 
which  the  use  of  NN  techniques  might  not  be  advisable.  Section  III  then  considers  the 
correlation  of  bearing  data  with  tracks  in  three-dimensional  space.  Section  IV  shows  how 
to  use  bearing  data  derived  from  sensors  with  overlapping  fields  of  view  to  determine  the 
three-dimensional  positions  of  the  objects  within  the  common  domain  of  observation.  Here 
the  pivotal  innovation  is  the  reduction  of  the  problem  from  enumerating  the  intersections 
of  bearing  vectors  recorded  by  different  sensors  to  evaluating  the  equality  of  pairs  of  special 
angles  computed  from  the  bearing  data.  This  format  facilitates  the  use  of  near-neighbor 
searches  and  fast  sorting  algorithms  and  vastly  reduces  the  computer  time  required. 

The  reader  should  note  that  the  term  “track”  might  be  used  here  in  a  less  precise  sense 
than  the  designers  of  correlator-tracker  software  might  desire.  Define  the  term  “target 
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map”  to  be  the  statistical  representation  of  a  state  vector  at  a  given  time  for  a  possibie 
object,  as  determined  from  sensor  contact  reports  [1].  A  track  is  then  a  time  sequence  of 
target  maps  corresponding  to  the  same  object.  The  statistical  nature  of  the  target  map 
derives  from  errors  implicit  in  the  measurement  and  estimation  process.  To  determine 
near  neighbors  in  the  present  paper,  the  authors  use  the  mean  locations  (or  expectation 
values  of  the  respective  positions)  of  the  contact  reports  and  target  maps.  Further,  the 
term  “track”  stands  for  the  mean  location  of  the  target  map  at  a  given  time  t. 
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II.  Data  Correlation  in  Many-Body  Scenarios 

In  this  discussion,  a  data  set  will  contain  “data  elements,”  each  of  which  corresponds  to 
one  of  a  set  of  objects  being  observed.  Thus  each  data  element  can  be  a  contact  report  from 
one  of  an  array  of  sensors  covering  a  region  of  physical  space  or  a  target  map  produced  by  a 
correlator- tracker  from  such  data.  The  statement  that  two  data  elements  are  “correlated” 
means  that  each  data  element  describes  or  pertains  to  the  same  physical  object  with  a 
likelihood  that  exceeds  some  threshold  value.  The  statement  that  two  disjoint  sets  of  data 
elements  are  “related”  means  that  each  of  the  sets  contains  at  least  one  data  element  that 
is  correlated  to  a  data  element  of  the  other. 

The  problem  considered  here  is  stated  as  follows:  given  M  >  2  related,  disjoint 
sets  of  data  (denoted  Si,  S2,  ...,Sm),  find  the  elements  of  the  respective  sets  which  are 
correlated  with  each  other.  That  is,  find  all  sets  {a,,  bj,  ...,  cm}  where  a,  £  Si,  bj  £ 
S2,  . . .,  and  cm  £  Sm,  so  that  cq,  bj,  ...,  cm  are  all  correlated.  “Data  correlation”  or 
simply  “correlation”  is  the  process  of  identifying  correlated  data  elements.  In  practice,  the 
determination  that  two  elements  of  different  data  sets  are  correlated  can  be  approximate 
(probability  of  correlation  less  than  one),  and  estimating  the  degree  of  correlation  requires 
the  calculation  of  a  correlation  score.  As  indicated  above,  one  example  is  the  act  of  finding 
the  track  (constructed  from  previous  data)  which  is  correlated  with  a  new  sensor  contact 
report.  Upon  establishing  this  relationship,  one  can  extend  the  track  to  later  times,  using 
the  new  contact  report  to  refine  the  trank. 

The  relevant  scenarios  comprise  tens  to  hundreds  of  thousands  of  rapidly  moving  ob¬ 
jects.  We  further  assume  that  the  period  under  consideration  is  a  few  hours  at  most. 
Because  of  the  large  numbers,  high  speeds  and  short  time  scales  involved,  time  delays 
of  even  a  minute  in  data  processing  could  be  too  large  to  react  to  the  situation.  One 
thus  immediately  concludes  that  data  correlation  in  this  case  requires  highly  parallel  com¬ 
puters  with  large  amounts  of  memory,  so  the  requisite  data  processing  can  be  as  nearly 
instantaneous  as  possible.  Only  analysis  and  numerical  simulations  of  realistic  scenarios 
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containing  such  large  numbers  of  objects  can  provide  estimates  of  the  computing  needs  in 
this  situation. 

II.  1  Why  Near-Neighbor  Algorithms? 

Work  has  begun  on  the  development  of  highly  parallel  correlation  algorithms  under  the 
assumption  that  scalar  algorithms  and  hardware  would  not  be  adequate  for  the  scenarios 
mentioned  above.  The  parallel  environment  would  permit  groups  of  correlation  calculations 
to  be  performed  simultaneously.  Ideally  the  individual  data  elements  involved  in  the 
parallel  calculation  correspond  to  the  same  observation  time  or  nearly  the  same  observation 
time  t0.  This  is  where  near-neighbor  algorithms  offer  the  greatest  advantage  over  brute 
force  methods. 

A  “brute  force”  method  computes  a  correlation  score  for  all  possible  sets 
{a*,  bj,  . . . ,  cm}.  When  the  datasets  Si,  S2,  . .  .,Sm,  contain  the  same  number  of  elements 
Ni  =  N2  =  . . .  =  Nm  2  N,  the  brute  force  calculation  scales  as  N v/ ,  where  M  is  the 
number  of  data  sets  to  be  correlated  at  one  time.  This  paper  considers  problems  for  which 
M  =  2  or  3. 

A  “near-neighbor”  algorithm  is  a  method  for  storing  and  accessing  data  on  N  objects 
so  that  the  following  conditions  hold: 

(1)  Storage  and  access  of  data  depends  on  object  parameters:  for  example,  position,  ve¬ 
locity,  or  attributes.  In  a  dynamical  system,  time  can  either  be  explicit  or  implicit  in  the 
near-neighbor  data  structure. 

(2)  The  storage  and  access  depend  on  indexing  the  data  or  pointers  to  the  data  (a  “data 
structure”)  so  that  the  indices  themselves  contain  information  on  parameter  values;  and 

(3)  On  average,  scaling  of  the  access  time  with  N  is  better  than  “brute  force"’  methods  by 
a  factor  of  approximately  N/  log  N  or  better. 

Near-neighbor  algorithms  eliminate  an  appreciable  fraction  of  the  combinations 
{a<,  bj,  ...,  cm}  from  consideration  before  correlation  scores  are  computed.  Consider 
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the  case  M  =  2.  For  each  a<  €  Si,  one  accesses  only  the  data  elements  in  S2  which 
satisfy  definite  criteria  for  being  near  neighbors  of  a,.  The  resulting  subset  of  S2  is 
NN2i(t)  =  {6j|  bj  NN  of  a<}  C  S3  and  the  correlation  score  is  computed  only  for  the 
near  neighbors  NN2i(i).  Depending  on  the  particular  NN  algorithm,  the  cost  of  accessing 
near  neighbors  for  each  ^  G  Si  scales  as  either  N°  (constant)  or  In N .  So  the  cost  for  N 
such  calculations  (all  elements  of  Si)  scales  as  N  or  N In .V,  both  of  which  are  less  than 
brute  force  scaling.  For  N  =  10, 000  objects  and  .V/  =  2,  the  use  of  NN  algorithms  gives  a 
potential  improvement  of  three  to  four  orders  of  magnitude  over  the  brute  force  approach. 
For  M  =  3,  the  improvement  can  be  seven  to  eight  orders  of  magnitude. 

To  summarize  the  above  discussion,  a  near-neighbor  algorithm  involves  the  construc¬ 
tion  and  use  of  an  intermediate  data  structure  which  permits  the  rapid  identification  of 
likely  candidates  (usually  a  relatively  small  number)  for  correlation  with  a  given  data  el¬ 
ement  before  the  actual  correlation  score  calculation  is  done.  If  the  same  near-neighbor 
data  structure  is  used  for  all  data  accesses,  the  scaling  with  N,  the  number  of  elements  in 
a  data  set,  is  faster  by  several  orders  of  magnitude  than  is  a  brute  force  algorithm,  when 
N  is  large. 

Table  1  lists  near  neighbor  algorithms  that  the  authors  have  developed  or  tested.  The 
quantity  N  is  the  number  of  elements  of  set  Si,  all  of  which  reside  in  the  near-neighbor 
data  structure.  Column  2  shows  the  scaling  of  the  time  required  to  set  up  the  NN  data 
structure.  The  scaling  of  the  cost  for  accessing  near  neighbors  of  a  given  element  of  Si 
appears  in  the  third  column.  For  this  operation,  NN  algorithms  are  superior  to  brute 
force.  On  the  other  hand,  column  2  indicates  that  if  only  a  few  accesses  occur  before 
the  data  structure  must  be  set  up  again  or  reconstructed,  the  total  cost  of  using  a  near 
neighbor  algorithm  could  scale  more  poorly  with  N  than  brute  force.  The  NN  algorithm 
is  not  necessarily  optimal  in  every  case.  For  use  of  a  NN  algorithm  to  be  beneficial,  the 
cost  of  organizing  the  data  according  to  a  NN  data  structure  and  computing  scores  for  the 
resulting  limited  subsets  of  Si  ®  S2  must  be  less  than  the  total  cost  (over  the  entire  set 
S2)  of  comparing  all  elements  bj  of  S3  with  all  elements  a{  of  Si. 
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Table  I  —  Scaling  of  Candidate  Gating  Algorithms 


*  Coefficients  can  vary  considerably 
t  Monotonic  Lagrangian  Grid 

*  Bilinear  Decomposition 


II. 2.  When  Will  NN  Algorithms  Work  Best  for  Tracking  and  Correlation? 

To  see  how  NN  algorithms  can  increase  the  speed  of  data  correlation,  one  can  analyze 
the  following  equation  for  the  cost  in  computer  time  of  correlating  all  JVn  elements  of  set 
S2  with  all  Ni  elements  of  Si: 

C  =  C*  +  C,  .  2.2.1 1 

In  Eq.(2.2.1),  Cn  is  the  cost  of  finding  the  elements  b}  £  S2  which  are  near  neighbors 
of  the  respective  a.{  £  Si;  and  Cs  is  the  time  required  to  compute  scores  for  all  of  the 
NN  pairings  €  NN2i(i),  for  i  =  1,  2,  . . . ,  N\.  The  scaling  of  each  term  with  Xi 

and  1V2  provides  a  basis  for  initial  comparison  of  NN  algorithms  with  brute  force  methods. 
The  constant  coefficients  multiplying  the  scaling  functions  will  vary  with  the  commuter 
hardware  and  the  degree  to  which  the  software  is  optimum  for  that  hardware. 

In  the  context  of  correlation  and  tracking,  define  set  Si  as  the  set  Sr  of  contact 
reports,  corresponding  in  the  most  general  case  to  different  observation  times  {tt  ;  i  = 

I, 2,...,  iVr},  where  N\  =  Nr.  the  number  of  contact  reports.  The  set  S2  is  then  equal  to 
St,  the  set  of  N?  =  Nt  tracks.  Thus  for  each  contact  report  a,  £  Si  our  objective  is  to 
identify  all  tracks  in  S2  that  are  correlated  with  a;.  The  tracks  in  S2  all  correspond  -o 
the  same  time  which  is  different  from  the  observation  times  for  the  set  {a,}  in  the  most 
general  case.  The  computation  of  a  correlation  score  requires  that  a  track  and  a  report 
correspond  to  the  same  time  [2],  so  that,  in  theory,  the  set  of  tracks  must  be  extrapolated 
to  the  time  corresponding  to  each  contact  report  a,  £  Si.  The  score  can  decrease  on 
average  with  track-report  distance,  but  is  not,  in  general,  a  function  of  the  distance  alone. 

II.  2.1.  “Brute  Force”  Methods 

A  brute  force  method  treats  each  contact  report  separately  from  the  others  and  per¬ 
forms  calculations  (not  necessarily  scores)  for  all  possible  pairings  of  tracks  with  contact 
reports.  Several  brute  force  procedures  Eire  possible.  Below  we  assume  that  extrapolation 
of  each  track  to  the  time  of  each  report  is  implicit  in  any  brute  force  algorithm. 
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Procedure  1.  For  each  contact  report  ax,  compute  the  “distances”  Dt]  to  ail  tracks  b}. 
For  each  track  satisfying  Z?,;  <  R,  where  R  is  a  distance  chosen  to  define  near  neighbors, 
compute  the  correlation  score.  The  scaling  of  the  cost  equation  (2.2.1)  with  jVr  and  Nt  is 

C*3  —  4b  V  V.  4-  4b  V  (9  9  9) 

'-'1  —  i*r  .iVf  -r  -4.lc  .yr  .  - ) 

In  this  equation,  the  coefficients  .4  axe  constant,  and  the  superscript  b  signifies  “brute 
force.”  The  coefficient  .4^.  depends  on  the  average  number  of  tracks  found  to  be  near 
neighbors  of  each  contact  report  and  the  cost  of  computing  a  score.  Because  the  first  step 
of  Procedure  1  culled  out  the  tracks  that  were  poor  candidates  for  correlation,  the  second 
term  does  not  depend  on  Nt.  The  coefficients  are  important,  since  the  actual  correlation 
calculation  will  be  much  more  expensive  than  the  calculation  of  distances  in  the  first  step. 
Notice  that  the  distances  Dxj  are  straightforward  to  calculate  only  when  physical  positions 
are  used  to  define  “near  neighbors”.  If  other  attributes  are  to  be  considered,  the  testing 
might  be  on  intervals  of  attribute  values  rather  than  on  a  normalized  “distance”  function. 

A  possible  drawback  of  Procedure  1,  as  described  above,  is  that  the  value  of  R  is 
somewhat  arbitrary  and  in  fact  might  have  to  be  inefficiently  large  to  ensure  finding  all 
the  correlations  with  acceptable  scores.  The  reason  is  that  the  score  is  not  a  function  of 
R  alone,  additionally  depending  on  the  statistical  distributions  defining  the  track  state 
vector  and  the  contact  report.  This  means  that  some  tracks  with  unacceptable  correlation 
values  would  be  selected  m  the  near-neighbor-access  step  of  Procedure  1.  Alternatively 
one  could  do  a  simple  comparison  to  find  all  tracks  within  an  initial  radial  interval  (0,  Rq) 
of  contact  report  a,-  and  compute  scores  for  this  subset  NN2i(i).  Then  choose  a  radius 
Ri  >  Rq.  For  the  remainder  of  the  tracks  (the  set  S2  -  NN2i(i)),  find  the  tracks  in  the 
interval  (Rq,R\),  compute  their  correlation  scores,  and  eliminate  those  tracks  from  the 
set  of  candidates  for  the  next  radial  interval.  This  process  continues  to  larger  distances 
Ri  until  all  correlation  scores  for  the  tracks  in  the  interval  (Ri,Ri+i)  fall  below  a  selected 
threshhold.  The  size  of  the  subset  of  S2  to  be  searched  and  the  associated  cost  decrease 
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with  successive  radial  intervals.  In  addition  the  algorithm  should  identify  all  correlations 
with  a  score  above  the  threshhold.  Testing  on  realistic  scenarios  is  necessary  to  determine 
the  efficiency  of  Procedure  1  relative  to  the  procedures  below. 

Procedure  2-  For  each  contact  report  a*,  compute  the  full  correlation  score  for  all  tracks. 
Keep  only  those  cracks  with  scores  above  some  minimum  acceptable  value  or  threshold. 
The  scaling  of  the  cost  equation  (2.2.1)  with  Nr  and  Nt  is 

c}  =  .4  Nr  N,  .  (2.2.3) 

Here  no  preliminary  determination  of  near  neighbors  has  occurred  before  the  correlation 
score  calculation,  so  that  Eq.(2.2.3)  indicates  greater  cost  than  both  parts  of  the  previous 
procedure,  given  that  calculating  the  correlation  score  is  the  most  expensive  part  of  the 
algorithm. 

Procedure  3.  For  each  contact  report  ai,  compute  and  sort  the  distances  Di}  to  all  tracks, 
and  compute  correlation  scores  in  order  of  increasing  distance  until  the  score  falls  below  a 
minimum  acceptable  value.  The  cost  of  this  is  harder  to  quantify  because  the  number  of 
scores  calculated  depends  on  a  score  threshhold  instead  of  a  distance  threshhold.  A  useful 
equation  is 


C3b  =  A*nNr  Nt  +  A*  Nr  Nt  In  Nt  +  -4b  /  Nt  .Vr  .  (2.2.4) 

The  coefficient  A^c  depends  on  the  cost  of  computing  a  score.  The  quantity  /  is  the  average 
fraction  of  the  Nt  tracks  for  which  correlation  scores  must  be  calculated  on  a  “per  contact 
report”  basis.  If  /  is  sufficiently  small,  this  procedure  should  cost  less  than  Procedure  2. 
In  addition,  if  /  ~  iVt_a,  where  a  >  0,  then  the  scaling  of  the  score  calculation  (third 
term)  is  better  than  in  Eq.(2.2.3).  This  method  also  eliminates  the  potential  drawback  of 
Procedure  1  that  the  value  chosen  for  R  is  somewhat  arbitrary  and  in  fact  might  have  to 
be  chosen  undesirably  large  to  ensure  finding  all  the  correlations  with  acceptable  values. 


9 


Notice  that  the  first  two  terms  have  potentially  more  expensive  scaling  than  the  third 
but  that  the  coefficients  Aj,,  and  IniV,  could  be  much  smaller  than  A^./.  In  addition, 
the  remote  possibility  exists  that,  even  for  large  numbers  of  objects,  the  coefficients  A 
could  be  more  important  than  the  scaling  with  Nr  and  Nt. 

Procedure  4.  Here  we  mention  a  nonhierarchical  cluster  algorithm  [ij,  which  is  discussed 
in  a  companion  document  [3]  and  analyzed  fully  there.  This  algorithm  identifies  ar.d  uses 
clusters  containing  tracks  and  reports  which  can  only  be  associated  with  members  of  the 
cluster  in  which  they  reside.  Based  on  this  information,  the  computation  of  correlation 
scores  has  scaling  ranging  between  N 3/2  and  iV2 ,  depending  almost  entirely  on  the  scenario. 
We  refer  the  reader  to  the  above  reference  [3]  for  further  information. 

An  important  factor  in  the  values  of  the  coefficents  for  the  above  procedures  is  the 
availability  and  type  of  computer  hardware.  The  speed-up  through  parallelization  and  the 
parallel  implementation  itself  will  likely  vary  with  the  particular  calculation  (each  term 
in  the  above  equations).  This  could  change  the  scaling  with  N  in  the  cost  equations 
or  change  the  relative  importance  of  the  individual  terms.  In  addition,  for  particuk  - 
multiple-hypothesis-tracking  (MHT)  algorithms,  some  of  the  procedures  above  might  no* 
be  possible  or  acceptable.  Most  importantly  for  the  brute  force  algorithms,  the  processing 
of  each  contact  report  will  be  the  same  whether  the  contact  report  time  stamps  {t,}  are 
different  or  not.  This  could  result  in  an  advantage  over  near-neighbor  algorithms,  which 
are  considerably  more  efficient  when  many  of  the  contact  reports  have  approximately  the 
same  time  stamps. 

II.2.2  Near-Neighbor  (NN)  Algorithms 

As  indicated  earlier,  NN  Algorithms  are  best  suited  for  processing  groups  of  contact 
reports  in  a  parallel  or  vector  fashion.  The  constraints  on  the  operational  sensor  system 
itself  can  limit  the  utility  of  such  algorithms,  for  example,  when  only  a  few  objects  are 
being  tracked  or  when  the  sensor  produces  reports  corresponding  to  different  times.  At 
least  three  cases  occur  in  correlator-tracker  problems: 
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Case  1.  Appreciable  numbers  of  contact  reports  in  Sr.  correspond  to  approximately  the 
same  observation  time  f0.  In  this  case,  the  same  near-neighbor  data  structure  containing 
the  tracks  extrapolated  to  tj  =  t0  can  be  used  to  process  those  contact  reports.  Then  the 
cost  equation  is 


cr  =  A"n  Nt  In  .V,  +  A?c  Nr  (2.2.5) 

or 

(W  -  (A?.y  JVf  +  (A?e)' iVP  ,  (2.2.6) 

depending  on  the  scaling  of  the  cost  for  setting  up  the  NN  data  structure.  Notice  that  in 
either  case,  the  scaling  with  Nr  and  Nt  is  superior  to  the  brute  force  methods.  Case  1 
applies  best  when  the  range  of  contact  report  observation  times  is  comparable  to  or  smaller 
than  the  minimum  time  for  any  track  to  overtake  another,  thus  disturbing  the  NN  data 
structure. 

Table  2  shows  timings  for  some  of  the  data  structures  in  Table  1  on  a  SUN  260  (scalar  ) 
workstation  for  the  case  that  all  of  the  contact  reports  have  the  same  time  stamp.  These 
tests  distributed  8K  sizeless  points,  designated  “tracks,”  and  8K  additional  sizeless  points, 
designated  “reports,”  over  a  three- dimensioned  space  uniformly  and  randomly.  The  “fc- d 
tree”  is  a  standard  ^-dimensioned  binary  tree  structure  and  search;  in  this  case,  k  =  3. 
The  various  methods  then  found  the  five  nearest  tracks  for  each  report.  We  see  that  near- 
neighbor  algorithms  provide  a  tremendous  advantage  over  brute  force  algorithms  when 
Case  1  applies. 

Case  2.  The  other  end  of  the  spectrum  is  the  situation  in  which  the  reports  correspond 
to  a  wide  range  of  observation  times,  all  widely  separated  from  each  other.  To  determine 
likely  track-report  pairings,  the  set  of  tracks  might  have  to  be  extrapolated  to  each  different 
contact  report  time  tj.  The  NN  data  structure  containing  the  tracks  might  then  have  to 
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Table  2  —  Preliminary  Tests  of  Candidate  Gating  Algorithms 


COSTS  MEASURED  IN  CENTRAL  PROCESSOR  SECONDS 
BILINEAR  DECOMPOSITION  (NEWLY  DEVELOPED  AT  NRL) 


be  set  up  anew  for  each  report.  Given  that  this  happens  approximately  Nr  times,  the  cost 
is 


C2n  =  A?„  Nr  Nt  In  Nt  +  A?e  Nr  (2.2.7) 

or 

(C?)'  =  (A?n)'  Nr  Nt  +  (A?e)#  -Vr  .  (2.2.3) 

Asstiming  that  Nr  ~  iVt  ~  jV,  the  actual  number  of  objects,  comparison  with  the  earlier 
equations  for  brute  force  methods  shows  that  in  this  case  the  scaling  for  NN  algorithms 
is  no  better  or  even  worse  than  brute  force.  Thus  if  the  contact  reports  are  processed  one 
at  a  time  and  the  observation  times  are  all  sufficiently  different,  NN  algorithms  might  not 
provide  an  advantage  over  brute  force  algorithms. 

Case  3.  An  intermediate  situation  prevails  if  the  various  objects  change  their  indices  in  the 
NN  data  structure  slowly  with  time.  This  means  that  their  relative  positions  in  parameter 
space  remain  the  same  over  appreciable  periods  of  time.  The  NN  data  structure  then 
remains  the  same  for  a  range  of  contact  report  observation  times  and  accessing  the  near¬ 
neighbors  of  a  contact  report  can  proceed  as  in  Case  1.  Of  course,  the  actual  spatial 
coordinates  associated  with  a  given  set  of  indices  in  the  NN  data  structure  do  change  and 
one  must  account  for  this  change  in  values  when  accessing  near  neighbors.  Methods  of 
doing  this  are  presently  under  development  by  the  authors  (e.g.,  [3]). 

Note  that,  as  in  Brute  Force  Procedure  1,  the  distance  R  or  the  parametric  intervals 
used  to  define  the  concept  of  “nearness”  might  have  to  be  larger  than  optimum  to  assure 
the  user  of  finding  all  candidate  pairs  with  correlation  scores  above  some  threshhold.  The 
choice  of  R  might  be  somewhat  arbitrary,  and  because  the  problem  is  inherently  statistical, 
the  chance  of  missing  a  required  candidate  on  the  first  pass  is  still  nonzero.  This  depends 
on  the  form  of  the  scoring  function.  As  in  brute  force  Procedure  1,  making  a  sequence  of 
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near-neighbor  search  and  scoring  procedures  at  succesively  larger  radius  might  overcome 
the  limitation: 

(1)  Using  a  near- neighbors  data  structure,  find  the  tracks  with  Dij  <  i?0>  for  some  rea¬ 
sonable  Rq  and  compute  scores.  Record  all  of  the  tracks  (subset  NN2i(i))  found  in  this 
step. 

(2)  Using  the  same  data  structure  and  a  larger  radius  Rx,  find  the  tracks  with  Rq  <  DXJ  < 
Ri  in  the  reduced  set  S2  -  NN2i(t).  Compute  scores  and  compare  to  the  threshhold. 
Eliminate  the  tracks  found  in  this  step  from  the  subset  of  tracks  subjected  to  subsequent 
searches. 

(3)  Repeat  step  (2)  in  the  interval  (i2x,f?2]-  Continue  this  procedure  for  subsequent  con¬ 
tiguous  intervals  until  all  scores  in  an  interval  fall  below  the  threshhold. 
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III.  Single  Sensor  Bearing  Data 

The  discussion  below  considers  cases  in  which  three-dimensional  tracks  are  to  be  cor¬ 
related  with  the  two-dimensional  bearing  data.  Figure  1  depicts  this  situation.  The  black 
dots  are  the  centra  of  the  various  “target  maps”  at  the  present  time.  Target  maps  repre¬ 
sent  (statistically)  the  state  vectors  of  potentially  real  objects  at  a  given  time  t.  Tracks  are 
then  time  sequences  of  target  maps,  each  sequence  corresponding  to  a  presumed  physical 
object.  Because  tracks  are  derived  from  incomplete  data  and  are  subject  to  ambiguous 
interpretations,  more  tracks  will  usually  exist  than  actual  objects.  As  indicated  earlier,  we 
sometimes  use  the  terms  “track”  and  “target  map”  interchangeably  for  convenience. 

The  coordinate  axes  at  the  left  side  of  Fig.  1  are  symbolic  of  any  universal  coordinate 
system  in  which  the  tracks  are  defined.  The  vector  S  defines  the  location  of  the  sensor  in 
this  coordinate  system.  Because  the  range  of  the  contact  reports  is  assumed  to  be  unknown, 
bearing  data  from  a  single  sensor  could  be  correlated  with  any  targets  lying  within  a  cone, 
whose  vertex  lies  at  the  position  of  the  sensor  and  whose  solid  angle  corresponds  to  the 
angular  resolution  of  the  sensor.  In  the  figure,  (A0,,A<£s)  define  the  allowed  angular 
interval  for  correlation  of  tracks  with  a  contact  report  at  the  spherical  coordinates  (04,  o,). 

In  this  section,  the  number  of  dimensions  of  the  parameter  space  defining  near  neigh¬ 
bors  is  n  =  2,  corresponding  to  {#,  <p),  and  the  data  arrays  in  computer  memory  will  be 
two-dimensional.  Picone  et  al.  [3]  discuss  examples  of  such  schemes,  two  of  which  are 
described  in  the  following  subsections. 

III.l.  Hockney’s  Method 

In  Hockney’s  method  [4],  an  n-dimensional  rectangular  grid  overlays  the  physical  space 
containing  the  tracks  and  contact  reports.  Each  side  of  a  rectangular  cell  corresponds  to 
one  of  the  physical  parameters  defining  the  correlation  calculation  (e.g.,  spatial  coordinates, 
velocity  components,  or  physical  attributes  of  the  objects  being  tracked)  and  has  length 
L\.  In  the  present  case,  l  is  one  of  the  angular  spherical  coordinates  centered  on  the  sensor 
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Fig.  1  —  Correlation  of  bearing  data,  represented  by  the  cone  emanating  from  the  sensor,  with 
tracks,  represented  by  black  circles.  Our  method  requires  the  transformation  of  the 
track  locations  to  a  coordinate  system  whose  origin  coincides  with  the  sensor  location. 


(9  or  <f>)  and  n  =  2.  To  find  which  rectangle  contains  a  given  point,  (0;,  4>i)  requires  two 

multiply  operations 

(/i.josckxi-l.Lx-p-l),  (3.i) 

L$j  l  Lj, 

where  (7),  Jj)  are  the  indices  of  the  cell  in  which  object  i  is  located  and  the  brackets  (jaj) 
indicate  the  greatest  integer  less  than  or  equal  to  the  quantity  a.  To  handle  cells  containing 
more  than  one  object,  one  employs  linked  lists.  Thus  the  near  neighbors  of  a  given  object 
are  those  within  the  same  cell  and  the  adjacent  cells.  The  cost  of  partitioning  the  objects 
among  the  cells  according  to  Eq.(3.1)  scales  as  Nt,  the  total  number  of  tracks.  This  method 
permits  parallel  calculations  of  correlation  scores,  but  should  not  permit  vectorization  on  a 
Cray  computer.  Because  the  indexing  of  the  objects  in  memory  bears  a  direct  relationship 
to  the  actual  (not  relative)  values  of  the  parameters  (0,  <j>)  for  each  object,  this  method  will 
permit  the  user  to  find  all  of  the  near  neighbors  corresponding  to  a  given  correlation  radius 
Rc.  For  nonuniform  distributions  of  objects,  many  cells  could  be  empty.  In  addition,  the 
values  of  Lg  and  are  somewhat  arbitrary.  Both  factors  could  cause  wasted  memory, 
and  the  latter  could  lead  to  inefficiency  and  total  costs  for  correlation  which  scale  more 
expensively  than  Eq.(2.2.6)  would  indicate.  In  parallel  processing,  load  balancing  becomes 
a  concern  because  some  nodes  can  have  many  neighbors  while  others  may  have  none. 

III.2.  The  Monotonic  Lagrangian  Grid  [3,  5] 

The  Monotonic  Lagrangian  Grid  (MLG)*  correlator  could  be  based  on  (0,  <j>)  in  a 
spherical  coordinate  system  having  its  origin  (r  =  0)  at  the  position  of  the  sensor.  The 
MLG  would  then  be  a  two-dimensional  data  structure  containing  both  the  bearings  of 
the  target  maps  relative  to  the  sensor  and  the  new  bearing  contact  reports.  The  pairs 
would  have  a  “monotonic”  order  in  the  sense  that 

*  Early  papers  used  the  term  “Monotonic  Logical  Grid”  for  this  data  structure. 
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9(i  +  l,j) 

<t>(i,j)<<t>(i,j  +  l)  ■ 

Then,  according  to  the  simplest  concept,  the  identification  of  potential  correlations  would 
only  require  computing  scores  of  target  maps  lying  within  index  neighborhoods  of  the 
contact  reports.  The  ordered  partition  of  Kim  and  Park  [6]  provides  a  way  to  augment  an 
MLG  so  that  a  tree  search  can  be  performed. 

The  calculation  of  the  MLG  indices  would  proceed  precisely  as  in  the  Cartesian  case  [3]. 
Transformation  of  the  target  maps  from  some  universal  coordinate  system  to  the  (r,  9,  cf>) 
sensor  system  would  require  only  highly  parallel  and  vectorizable  calculations  which  would 
also  be  necessary  for  any  non-MLG  technique!  This  process  scales  as  order  ( Nt )  or  order 
(1),  depending  on  the  computer  architecture  and  memory  and  the  number  of  target  maps 
to  be  processed.  For  scalar,  MIMD  (multiple  instruction  stream,  multiple  data  stream)  or 
pipeline  architectures,  the  MLG  bearing  data  structure  might  consist  of  only  pointers  to 
the  3D  MLG  target  map  data  array/file  and  the  sensor  contact  report  array/file.  Thus  the 
contact  reports  and  target  maps  need  not  be  stored  twice.  For  SIMD  (single  instruction 
stream)  architectures,  actually  combining  and  storing  the  data  in  a  contiguous  region  of 
memory  might  be  optimal. 

Because  the  MLG  is  determined  primarily  by  the  relative  locations  of  the  target  maps 
and  not  the  origin  of  the  coordinate  system,  the  motion  of  the  sensor  should  not  have  much 
of  an  effect  on  the  structure  of  the  MLG.  Only  local  swaps  of  the  target  maps  in  the  data 
structure  would  be  necessary  to  maintain  the  MLG  over  the  time  interval  during  which 
new  sensor  data  are  processed,  if  this  is  at  all  necessary.  Our  studies  have  shown  such 
restructuring  to  be  very  fast.  In  passing,  we  note  that  the  motion  of  the  sensor  relative 
to  another  (universal)  coordinate  system  would  not  affect  the  speed  of  transformation  of 
spatial  coordinates  from  one  system  to  the  other.  This  is  so  because  the  motion  merely 
changes  the  value  of  the  vector  S  relating  the  origins  of  the  two  systems. 
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A  drawback  to  the  MLG  for  CT  problems  is  that  the  storage  of  data  in  computer 
memory  can  result  in  local  distortions  of  the  relative  configuration  of  objects  in  physical 
space.  If  the  same  index  offsets  are  used  to  access  the  tracks  that  axe  near  neighbors  of 
each  contact  report,  some  of  the  near  neighbors  can  be  missed  if  these  index  offsets  axe 
chosen  to  be  too  small.  Thus  large  index  offsets  might  be  necessary  in  dense  regions  to 
insure  that  all  of  the  desired  neighbors  are  found.  This  would  mean  that  many  undesirable 
candidates  would  be  found  as  well,  increasing  the  number  of  correlation  calculations  above 
the  optimum  level.  This  reflects  a  corresponding  situation  of  many  contact  reports  in 
a  single  cell  of  the  Hockney  data  structure  or  many  empty  cells,  if  the  cells  are  small. 
Apparently  some  information  about  the  actual  parameter  values  corresponding  to  MLG 
indices  would  be  necessary  to  improve  MLG  performance  under  such  conditions.  This  is 
the  essence  of  the  fc-nearest  neighbor  finding  algorithm  of  Kim  and  Park  [6],  which  can 
combined  with  an  MLG  data  structure. 

We  assume  that  the  reports  all  correspond  to  approximately  the  same  time  t0.  If  this 
is  not  the  case,  one  might  still  be  able  to  use  the  procedure  given  below  with  adjusted 
values  of  Ad  and  A <f>,  the  angular  uncertainties  in  the  observations  (step  3).  An  outline  of 
the  MLG  algorithm  is  as  follows: 

(1)  If  necessary,  compute  the  positions  of  the  target  maps  (T.)  relative  to  the  location 
of  sensor  S  and  corresponding  to  the  time  t0  at  which  the  observations  took  place. 
Derive  for  each  map  T;.  Referring  to  Fig.  1,  one  obtains  the  position  of  track 

i: 


where 


RJf  =s  R,  -  S, 


(3.2.1) 


Ri.  ->  ( n,di ,  (pi). 


(3.2.2) 
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(2)  Define  a  convenient  2D  MLG  containing  all  (JVP)  contact  reports  received  from  sensor 
5  and  the  angular  coordinates  of  all  (iVt)  candidate  target  maps  at  the  corresponding 
time  t0 •  Here  the  subscript  “f”  signifies  “tracks.”  Ideally,  the  product  of  the  dimensions 
N,  and  Nj  of  the  MLG  would  equal  the  total  number  of  angle  pairs  .Vr  4-  .V,  =  .Vr(. 
However,  .Vrf  might  not  be  factorizabie  into  che  MLG  dimensions  desired  see  text 
below).  In  Eq.(3.2.3)  below,  we,  therefore,  allow  for  the  addition  of  A.V  "hoies" 
to  make  up  this  difference.  A  hole  is  an  element  of  the  MLG  which  has  specific, 
convenient  values  of  (0fc,<pii)»  but  does  not  correspond  to  a  track  or  a  report.  The 
MLG  data  structure  thus  contains  the  set 

Ml  =  {(■£  jki  yjfc)|  (xjki  Ujk  )£{(^»  t  d>)  }  hi  {(#a, ©j  J }  U  ^  ^  ~($h ,  Oh  j; 

h  td.-.d) 

xjk  ^  x(j+l)ki  Vjk  5:  Vj(k  + 1)}  • 

(3)  Search  appropriate  index  neighborhoods  of  sensor  bearing  contact  reports  {(0*,  Oj)} 
for  sufficiently  correlated  target  maps.  The  degree  of  correlation  may  be  defined  by 
the  requirements  that 


\6.-6i\  <  A 9rt 
| d«  —  di  |  <  Adrt 


(3.2.4) 


where  A0r<  and  Adrt  are  measures  of  the  angular  resolution  of  the  sensor  or  a  similar  factor 
defining  the  required  degree  of  agreement  between  contact  reports  and  track  bearings. 


Item  (2)  above  includes  the  possibility  that  M2  would  require  a  few  more  (A.V) 
elements  than  Nrt  in  order  that  the  total  number  of  elements  N\i  is  factorizabie  in  the 
desired  manner.  For  example,  suppose  we  define  the  two  dimensions  of  M2  by 


Nj  =  Nk  =  [VN7t]  +  1, 


(3.2.5) ' 


where  [1]  signifies  the  greatest  integer  less  than  x.  Then  the  quantity 
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(3.2.6) 


A  N  =  N]  -  Nrt 

is  the  extra  number  of  empty  memory  locations  or  “holes”  within  M2 ,  In  this  case  the  holes 
are  just  a  convenience  in  setting  up  the  MLG  and  need  not  participate  in  the  correlation 
calculation.  Choosing  hole  locations  in  regions  with  a  statistical  paucity  of  real  data  might 
appreciably  regularize  the  MLG  structure. 
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IV.  Three  Sensors  with  Overlapping  Fields  of  View 

When  one  has  multiple  sensors  with  overlapping  fields  of  view,  the  bearing  data  can 
provide  enough  information  to  derive  the  positions  of  the  objects  being  observed  without 
reference  to  tracks.  The  larger  the  number  of  sensors  involved,  the  fewer  will  be  the 
pathological  cases  for  which  ambiguous  results  are  obtained.  In  the  brute  force  approach, 
the  bearing  data  may  be  represented  graphically  as  rays  (or  cones)  radiating  from  each 
sensor,  and  the  intersections  represent  possible  positions  of  objects  being  observed.  For 
m,  sensors  and  Nr  contact  reports  per  sensor,  checking  the  contact  reports  sequentially  for 
intersections  is  very  expensive  from  a  combinatorial  point  of  view,  scaling  approximately 
as  N™’ . 

A  method  for  finding  these  intersections  (correlations)  based  on  neax-neighbor  or 
sorting  techniques,  however,  scales  as  iVrlniVr,  which  is  a  major  improvement  over  brute 
force  for  m,  =  3  or  higher.  Using  highly  parallel  computer  architectures  could  cut  this 
scaling  by  a  factor  of  Nr.  This  section  describes  such  an  approach  for  the  problem  of  three 
sensors,  each  of  which  provides  bearing  data  only.  Our  interest  is  in  the  objects  being 
observed  in  a  region  of  space  viewed  by  all  three  sensors.  (As  this  report  was  going  to 
press,  the  authors  received  information  indicating  that  the  algorithm  presented  below  was 
developed  earlier  by  R.  Varad  and  perhaps  others  [7].  We  requested  reprints  of  the  papers 
in  question  but  have  not  yet  received  them  to  confirm  the  prior  claim.) 

Figure  2  shows  a  schematic  diagram  of  this  problem.  The  three  sensors,  labelled 
by  A,  B ,  and  C,  produce  sets  of  contact  reports  that  are  subscripted  respectively  by 
a,  6,  and  c,  and  labelled  respectively  by  i,  j,  and  k.  Originally  the  reports  give  the 
spherical  angles  (9,  <j>)  defining  a  bearing  vector,  denoted  by  eaij,  relative  to  a  specific 
sensor.  The  subscript  a  gives  the  sensor  indentification  while  p  labels  the  particular 
contact  report  produced  by  sensor  a.  Equation  (4.1)  represents  this  progression: 
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Fig.  2 


Conversion  of  bearings,  produced  by  multiple  passive  sensors  with  a  common  field  of 
view,  to  three-dimensional  positions.  The  sensors  are  labelled  .4,  jE?,  and  C .  and  the 
objects  being  observed  are  represented  by  black  circles. 
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(4.1) 


A  : 

eai  —  <Pa(*))i  *  —  h  2,... 

,  Na 

B  : 

ebj  =  (9b(j),<t>bU))i  j  =  1.  2,.. 

N„ 

C  : 

ec*  =  (0c(A:),<Mfc));  k  =  l,  2,. 

Ne 

For  the  purpose  of  discussion,  assume  that  the  sensors  produce  exact  bearing  data.  Then, 
as  Fig.  2  shows,  an  object  observed  by  ail  three  sensors  will  reside  at  the  intersections  of 
the  bearing  vectors  eai,  e&y,  and  ecjt,  for  specific  values  of  i,  j,  and  k.  Notice  that  the 
indices  i,  j,  and  k,  corresponding  to  a  threeway  intersection,  also  provide  unique  labels 
for  the  physical  object.  As  indicated  above,  the  number  of  possible  intersections  is  the 
product  Na  x  JVj  x  Nc  or  TV3,  if  Na  =  N(,  =  Nc  =  N.  Because  this  scaling  produces  a 
“combinatorial  explosion”  as  the  number  of  contact  reports  increases,  a  better  method  of 
correlating  the  sensor  data  is  necessary  to  derive  positions  for  the  objects  being  observed. 
The  following  discussion  describes  such  a  method. 

TV. 1.  Transformation  of  Bearing  Data  to  New  Angles 

Our  solution  of  this  problem  is  to  convert  the  spherical  coordinates  {9,  <f>)  to  two  new 
angles  which  define  the  planes  in  which  the  bearing  vectors  {ea,}  lie  relative  to  the  plane 
containing  the  three  sensors.  A  further  assumption  then  is  that  the  locations  of  the  three 
sensors  are  also  known. 

For  the  purpose  of  introducing  the  method,  our  discussion  does  not  include  the  angular 
resolution  of  the  sensors.  Errors  in  the  angular  data  produced  by  the  sensors  translate  into 
errors  in  the  new  angles  resulting  from  the  transformation.  The  Appendix  presents  a  prac¬ 
tical  implementation  of  the  ideas  presented  here  and  discusses  effects  of  sensor  resolution. 
The  next  subsection  on  practical  algorithms  and  tests  includes  a  simple  representation  of 
the  angular  error.  Another  ingredient  of  errors  in  the  angles  is  differences  in  observation 
times  of  the  contact  reports  produced  by  the  sensors.  Here  we  do  not  compute  the  effects 
of  differences  in  “time  stamps”  of  the  reports. 
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Fig.  3—  Conversion  of  the  bearing  angles  (Ob(j),  <i>b{]))  =  e&;,  corresponding  to  observation 
j  by  sensor  B,  to  new  angles  $*c(j)  and  The  new  angles  correspond  to  the 

planes  defined  by  the  combination  of  e^  with  each  of  the  lines  between  sensor  B  and 
the  other  sensors.  Notice  that  the  lines  of  sight  for  defining  the  angles  are  from  sensor 
B  toward  sensor  C  and  from  A  to  B ,  according  to  the  convention  in  Section  IV.  1  of 
the  text.  The  sensor  whose  letter  is  in  parentheses  is  hidden  from  the  observer  by  the 
other  sensor. 
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Now  consider  the  bearing  vector  ejy  for  report  j  of  sensor  B.  This  situation  is  repro¬ 
duced  in  Fig.  3a,  without  the  other  objects  and  bearing  vectors  in  Fig.  2.  Notice  in  Fig. 
3a  that  the  bearing  vector  ejj  in  combination  with  the  vector  e&c  along  the  line  joining 
sensors  B  and  C  defines  a  unique  plane.  Similarly  in  combination  with  ea6  defines 
another  unique  plane.  Further,  these  planes  intersect  at  the  bearing  vector  ej,;,  so  that 
knowledge  of  the  planes  is  equivalent  to  the  knowledge  of  the  bearing  vector. 

The  next  step  is  to  define  each  of  the  above  planes  in  terms  of  a  unique  angular 
coordinate.  Figures  3b  and  3c  indicate  how  this  is  done.  Figure  3b  shows  the  situation  as 
viewed  by  an  observer  whose  line  of  sight  lies  in  the  plane  of  the  sensors  and  is  directed 
along  the  line  connecting  sensors  B  and  C  (vector  eic).  The  plane  of  the  sensors  appears  as 
a  horizontal  line  to  the  observer.  The  line  connecting  B  and  C  appears  to  the  observer  as 
a  point  at  the  location  of  sensor  B.  The  plane  BCj  defined  by  the  bearing  vector  ejj  and 
the  vector  e&c  appears  as  a  line  intersecting  the  plane  of  the  sensors.  The  angle  of  rotation 
*6,0')  between  the  plane  BCj  and  the  plane  ABC  containing  the  sensors  is  unique  and 
thus  defines  the  plane  BCj.  Now  view  the  same  system  along  the  line  from  sensor  A  to 
sensor  B ;  the  observer  sees  the  situation  depicted  in  Fig.  3c.  In  the  same  manner,  the 
angle  ^b(j)  uniquely  defines  the  plane  ABj. 

This  method  has  thus  converted  the  angles  (&i(j),  4>b(j))  =  to  two  new  angles 
(^4c0)>  ^ai(i))*  The  new  angles  correspond  to  the  planes  defined  by  the  bearing  vector 
ebj  and  the  lines  joining  sensor  B  with  the  other  two  sensors.  These  planes  intersect  to 
form  the  bearing  vector  ebj.  The  same  procedure  performed  on  all  bearing  vectors  from 
sensor  A  provides  the  set  {($“4(i),  ^co(O)!  *  =  1>  2,...  JVB },  and  similarly,  the  bearing 
vectors  from  sensor  C  provide  the  set  {($4c(k),  $*„(&));  i  =  1,  2, . . .  Nc}.  By  convention, 
the  observer  views  the  system  in  the  plane  and  parallel  to  the  lines  joining  pairs  of  sensors 
as  follows: 

(1)  For  sensors  A  and  B,  along  the  direction  from  .4  to  B ; 

(2)  For  sensors  B  and  C,  along  the  direction  from  B  to  C;  and 
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(3)  For  sensors  A  and  C,  along  the  direction  from  C  to  A. 

This  maintains  a  cyclic  order  among  the  triplet  A  —  B  —  C.  The  index  i,  j,  or  k  and  the 
bold  superscript  a,  b,  or  c,  respectively,  tells  us  which  sensor  has  produced  a  given  angle 
$  in  the  above  sets. 

Now  return  to  Fig.  2  and  the  case  in  which  all  three  sensors  have  observed  the  same 
object.  Let  us  assume  that  the  bearing  vectors  eai,  e4j,  and  ec^  intersect  at  the  location 
of  an  object  as  shown.  First  consider  eat  and  e4j.  Because  the  two  vectors  intersect,  they 
lie  in  the  same  plane.  Further  the  line  AB  also  lies  in  the  same  plane.  Thus  the  plane 
ABi  corresponding  to  angle  <$“4(i)  (derived  from  sensor  A)  is  equivalent  to  the  plane  ABj 
corresponding  to  angle  (derived  from  sensor  B).  The  condition  that  the  bearing 

vectors  eai  and  e4j  correspond  to  the  same  object  (i.e.,  that  the  vectors  intersect)  is  thus 


*:.(•')- «rtO)  (4.1.1a) 

Similar  reasoning  for  the  bearing  vectors  from  the  other  pairs  of  sensors  B  —  C  and 
A  —  C  gives  us  the  conditions: 


(4.1.16) 

(4.i.ia) 

Summarizing  the  above  is  the  following  rule  for  determining  all  triplets  of  indices 
(i,j,k)  which  correpond  to  bearings  from  sensors  .4,  B,  and  C  that  mutually  intersect: 

Rule:  The  bearings  from  sensors  A ,  B,  and  C  with  respective  indices  i,  j,  and  k 
mutually  intersect  at  the  same  point  if  and  only  if  Eqs. (4.1.1)  are  satisfied. 

Thus  one  can  find  the  positions  of  the  objects  from  bearings  produced  by  three  sensors 
with  a  common  field  of  view  by  the  following  steps: 
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(1)  For  each  sensor,  transform  the  spherical  coordinates  for  each  data  point  to  the  angles 
<$,  as  defined  above. 

(2)  Search  the  sets  of  transformed  angles,  for  those  which  satisfy  Eqs.  (4.1.1). 

(3)  Compute  the  intersection  location  for  each  triplet  of  bearings  found  in  step  2. 

This  formulation  permits  the  use  of  sorting  or  near-neighbor  algorithms  to  find  the 
( i,j,k )  corresponding  to  the  same  object.  For  example,  to  find  those  i  and  j  satisfying 
Eq.(4.1.1a),  combine  the  set  {$“*(0}  coming  from  sensor  .4  with  the  set  {$^bU)}  from 
sensor  B  and  then  sort  the  entire  set  according  to  the  values  of  these  angles.  Those  angles 
that  are  equal  will  be  near  neighbors  in  the  sorted  data  array  and  finding  all  such  pairs 
will  scale  as  N,  the  number  of  elements  in  the  data  structure.  The  next  section  describes 
testa  on  various  sorting  and  near-neighbor  methods  designed  to  find  the  data  satisfying 
Eqs. (4.1.1). 

IV.2.  Representative  Algorithms  and  Tests 

The  algorithms  and  tests  below  account  for  errors  in  the  bearing  data  due  to  finite  sen¬ 
sor  resolution.  Using  the  notation  of  the  previous  subsection,  the  three  sensors  A,  B,  and  C 
each  produce  bearing  data  transformed  to  the  respective  sets  of  angle  pairs 

Sensor  A  :  Sa  =  {($“6(i),  $“<.(0)  =  ea,  |  «  =  1,  2, . . . ,  Na} 

Sensor  B  :  Sb  =  {($£k(j),  $ic(j))  =  ekj  |  j  =  1,  2, . . . ,  Nb}  (4.2.1) 

Sensor  C  :  Sc  =  {($?<.(£),  $“,(*))  =  eck  \  k  =  1,  2, . . . ,  Nc} 

Given  Sa,  Sb,  and  Sc  as  input,  we  will  examine  three  algorithms  for  finding  A  = 
{(i,  j.  A:)  |  such  that  e„j,  ejy,  ecfc  satisfy  Eq.(4.1.1)}.  Thus  A  is  the  set  of  index  triplets 
that  correspond  to  bearing  vectors  (one  from  each  sensor)  which  intersect  at  the  same 
“point”  (actually  a  region  defined  by  the  errors  of  the  sensors).  Given  the  bearing  vectors 
i,  e& j,  and  ec*,  one  can  directly  compute  the  coordinates  (x,  y,  z )  of  the  potential  object. 
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In  the  tests  below,  the  quantity  e  denotes  the  characteristic  error  in  the  bearing  data  and 
is  a  nonzero  parameter. 

IV. 2.1.  Parallel  Matching  Method 

This  method  requires  two  one-dimensional  sorts  of  the  angle  data  and  thus  scales  as 
N  In  N.  The  term  “parallel”  here  refers  to  the  fact  that  the  first  two  steps  of  the  procedure 
may  be  performed  simultaneously  on  a  minimum  of  two  processors.  Each  of  the  major  steps 
of  the  procedure  clearly  contains  similar  searches  whose  lengths  depends  on  the  data  sets. 
The  degree  of  parallellization  will  depend  on  the  architecture  of  the  particular  computer 
being  used.  The  steps  are  as  follows: 

(1)  Form  a  set,  denoted  Sab,  by  combining  the  angles  $“6(i)  6  Sa  with  the  angles  b(j )  € 
Sb-  Sort  Sab  and  search  outward  from  each  element  until  the  difference  between 

angle  values  exceeds  e.  This  gives  us  the  following  subset  of  Sab  which  contains  the  solution 
along  with  some  false  solutions: 


su  =  {su(o  =  u  i  *:»w  - « £  < «;»« + <}|*' 


(4.2.2) 


(2)  Perform  the  same  operation  with  the  angles  G  Sa  and  $*„(£)  6  Sc  in  a  set 

denoted  Sca.  This  gives  us  the  subset  of  Sca: 


su  =  {s;.(0  =  {(*  i  *;.(•)  -  s  <  *;.(»)  <  + <>|; 


1.2 .  .V.}  .  (4.2.3) 


(3)  Now  for  each  i,  compare  all  $4c(j)  for  j  G  Sab(i)  with  all  $4c(fc)  for  k  G  S*a(i).  This 
gives  us  the  index  triplets  ( i,j,k )  that  we  seek. 

These  subsets  Sab  and  should  contain  few  bearings  relative  to  the  total  number  pro¬ 
duced  by  a  sensor,  assuming  that  the  sensor  resolution  e  is  small  relative  to  the  average 

interobject  angular  separation.  In  fact,  we  have  assumed  that  the  subsets  Sab  and  S'a  are 
small  relative  to  the  binary  logarithms  of  the  number  of  elements  in  their  parent  sets  Sab 

and  Sca.  Otherwise  the  search  for  the  correlated  elements  making  up  set  Sab,  for  example, 
would  be  done  by  constructing  a  binary  tree  containing  all  $“(,(1)  €  Sa  and  then  searching 
the  tree  for  the  $“4(i)  correlated  with  each  e  Sb- 
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IV.2.2.  The  Monotonic  Lagrangian  Grid 

The  first  step  performed  by  the  MLG  method  is  to  construct  the  set  S‘ b  exactly  as  in 
the  parallel  matching  method.  For  each  i,  the  angles  $“6(i)  6  Sa  and  $kfc(j)  for  j  €  S^b(i) 
will  then  satisfy  Eq.(4.1.1a)  within  the  accuracy  defined  by  e.  Now  take  the  set  of  their 
counterpart  angles 

Sfb  =  {($:.(0,*»e0-))  I  i  =  1,  2,...,  Va;  je  s;b(*)}  (4.2.4) 

and  combine  these  pairs  with  their  counterpart  pairs  ($£„(&),  $£c(k))  €  Sc  produced  by 
Sensor  C  in  a  two-dimensional  MLG  data  structure: 

A  =  Sfb  U  Sc  .  (4.2.5) 

The  four  angles  satisfying  the  remaining  two  equations  (4.1.1b)  and  (4.1.1c)  should  be 
found  in  proximate  pairs  in  A  with  one  pair  from  each  of  the  constituent  sets  in  Eq.(4.2.5). 
As  noted  above,  when  one  searches  small  fixed-size  neighborhoods  of  each  point  in  the 
MLG,  the  search  is  rapid  but  some  of  the  correlations  might  be  missed.  Thus  one  first 
finds  a  high  percentage  of  the  solutions,  as  shown  by  the  data  in  Section  IV. 3.  Then 
one  can  repeat  the  procedure  on  a  much  smaller  data  set  consisting  of  the  remaining  few 
percent  of  the  contact  reports. 

IV.2.3.  A  Shortcut  Method 

The  value  of  this  shortcut  method  lies  in  finding  an  acceptable  percentage  of  possible 
correlations  rather  than  all  of  them  while  reducing  the  amount  of  computation  significantly. 
As  in  the  previous  methods,  compute  the  set  S^b  defined  by  Eq.(4.2.4).  For  each  i,  find 
the  element  j  €  S*b(z)  which  minimizes  the  deviation  of  $^t(ji)  from  $“6(i).  Then  vise  the 
triplet  ($*6(i),  $4C(j),  $*a  (i))  as  the  solution.  The  philosophy  behind  this  approach  is  that 
if  e  is  small,  then  S'b  will  also  be  small.  The  best  value  of  j  will  then  have  a  high  likelihood 
of  representing  true  correlations  even  without  confirmation  from  sensor  C.  Thus,  unlike  the 
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MLG  method  which  simply  fails  to  report  some  correlations,  the  inaccuracy  of  the  shortcut 
method  results  from  the  fact  that  some  false  reports  may  be  generated.  Another  variation 
of  this  method  might  be  to  report  every  triplet  ($* 6(i),  „(*'))  35  a  correlation  for 

t  =  1,  2, . . . ,  Na  and  j  6  S'b(i).  This  would  ensure  that  in  addition  to  some  false  reports, 
a  report  would  be  generated  for  each  of  the  real  objects  observed. 

IV.3.  Performance  Results  on  a  SUN-260  Workstation 

The  following  tables  provide  relative  performance  information  about  the  three  al¬ 
gorithms  described.  The  tests  took  place  on  a  SUN-260  Workstation,  which  is  a  scalar 
machine.  A  floating  point  accelerator  enhanced  the  speed.  Table  3  displays  execution 
times  and  accuracy  for  Na  =  JVj  =  JVC  =  N  €  {1 6K,  Z2K, . . . ,  128 AT}  objects,  with  each 
test  including  an  additional  10%  spurious  sensor  readings.  The  value  of  e  was  chosen  in 
each  test  so  that  A AB(i),  the  average  number  of  elements  in  Sab(i),  was  approximately 
1.1.  Table  4  presents  the  results  of  tests  performed  with  fixed  N  =  16 JU  and  varying  e 
(or  A This  table  reveals  how  the  size  of  Sab(i)  (as  well  as  S'a(i))  for  the  paral¬ 
lel  matching  method)  affects  the  performance  of  each  algorithm.  Figure  4  displays  the 
execution-time  breakdown  of  the  parallel  matching  and  MLG  methods  for  N  —  128A". 

The  data  for  all  tests  were  obtained  using  a  random  number  generator  rand()  which 
produces  a  uniform  distribution  of  values  between  0  (inclusive)  and  1.  The  first  N  values 
for  each  sensor  were  generated  as  follows  to  ensure  N  correlations: 

for  i  =  1  to  N  and  i  =  j  =  fc; 

$?6(0  =  rand();  ^b(j)  =  $*6(j’)  ±rand()  ■  e; 

$“a(0  =  rand();  $eca(k)  =  $‘a(i)  ±rand()  •  «; 

$^0)  =  rand();  $»c(£)  =  $£.(j)  ±rand()  •  e; 

end; 

A  random  number  of  spurious  readings  (amounting  to  around  10%  of  N)  were  then 
added  to  each  sensor  and  the  order  of  the  readings  permuted  so  that  eai,  ejj,  and  eck 
would  not  necessarily  correlate  for  i  =  j  =  k  .  For  simple  statistical  reasons,  none  of 
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Table  3  —  Timing  of  Bearing  Data  Fusion  Algorithm  vs.  /V 


RESULTS  FOR  n  =  16,000 


Correlation  Method  I  Time  (secs)  |  Correlations 


8.60 


11.33 


12.95 


MLG  (offset  =  3)  !  15.10  99.6% 


3.95 


Parallel  Matchin 


MLG  (offset  »  2) 


Shortcut  Method 


100.0% 


99.5% 


99.6% 


99.6% 


98.0% 


RESULTS  FOR  n  =  32.000 


Table  4  —  Timing  of  Bearing  Data  Fusion  Algorithm  vs.  Average  Data  Overlap 


|AA*  (/)  I 


Correlation  Method 


Parallel  Matching 


MLG  (of&et  =  1) 


Shortcut  Method 


=  1.1  (n  =  16000) 


Time  (secs) 

Correlations 

8.48 

100.0% 

11.15 

99.4% 

3.98 

98.0%  i 

|A  A“  (i )  I  =  1.2  (n  =  16000) 


Correlation  Method 

Time  (secs) 

Correlations 

Parallel  Matching 

8.79 

100.0% 

MLG  (of&et  *  1) 

11.34 

99.4% 

Shortcut  Method 

3.99 

95.8% 

|AAa  (i)l 


Correlation  Method 


Parallel  Matching 


=  1.8  (n  =  16000) 


Time  (secs)  Correlations 


9.33  100.0% 
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Execut ion-T  ime  Breakdown 
For  n  =  128,000 


Ml g  Method  Parallel  Matching  Method 

Fig.  4  —  Execution  time  breakdowns  for  the  MLG  and  the  Parallel  Matching  Method.  The 
notations  “abba,”  “acca,”  and  “bccb”  correspond  to  steps  1,  2,  and  3  of  the  latter 
method,  respectively. 
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the  spurious  readings  produced  false  correlations  in  the  tests  performed.  To  produce  false 
correlations  using  these  algorithms,  much  larger  values  of  sensor  error  t  or  many  more 
contact  reports  would  be  required. 

IV. 4.  Conclusions  from  Tests 

If  the  accuracy  of  the  shorcut  method  is  acceptable  for  a  particular  application,  'hen 
execution  speed  clearly  makes  it  the  method  of  choice;  however,  the  high  sensitivity  to 
measurement  error  (i.e.,  e)  as  target  densities  increase  probably  makes  the  shortcut  method 
impractical  for  use  in  a  tracking  system.  This  being  the  case,  the  parallel  matching  method 
must  be  considered  the  superior  algorithm  on  the  SUN-260  because  it  is  100%  accurate 
and  faster  than  the  MLG  method.  Relative  to  brute  force  algorithms,  either  method  is 
much  faster  because  of  superior  scaling  with  the  size  of  the  data  set. 

TV. 5  Estimates  of  Performance  on  Other  Hardware 

Table  5  gives  estimates  of  the  performance  of  our  multisensor  fusion  algorithm  on 
various  machines,  both  sequential  and  parallel.  The  numbers  are  based  on  tests  using  the 
MLG  and  on  relative  speeds  of  the  computers.  Under  each  machine  name  is  the  source 
of  the  estimate.  The  timing  estimates  shown  correspond  to  32.000  objects  as  viewed  by 
the  three  sensors.  On  various  vector  and  parallel  machines,  the  processing  time  is  on  the 
order  of  1  cpu  second.  The  algorithm  has  thus  reduced  the  complexity  of  the  problem  to 
the  point  at  which  existing  hardware  can  process  the  data  rapidly  for  large  numbers  of 
objects.  With  expected  improvements  in  processing  speed  over  the  next  few  years,  data 
processing  systems  should  easily  meet  the  requirements  for  bearing  to  position  conversion 
on  large  numbers  of  objects  with  realistic  values  of  sensor  false  alarm  and  error  rates. 
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Table  5  —  Estimate!)  of  Performance  of  Rearing  Data  Fusion  Algorithm  on  Various  liurilwaic 
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V.  Summary 

This  report  has  presented  a  preliminary  view  of  the  use  of  neax-neighbor  algorithms 
in  the  processing  of  data  on  the  bearings  of  the  objects  being  observed.  By  assumption, 
the  sensors,  such  as  infrared  imaging  devices,  have  provided  either  poor  range  information 
or  none  at  all.  Our  orientation  has  been  toward  correlation  and  tracking  of  large  numbers 
of  high-speed  objects.  The  purpose  of  the  near-neighbor  algorithms  has  been  to  organize 
the  data  on  a  set  of  objects  to  permit  rapid  access  of  the  data  with  a  high  correlation  to 
a  given  member  of  another  set.  The  combinatorially  inefficient  “brute  force'1  method,  in 
which  all  possible  correlations  between  the  two  sets  are  considered,  provides  a  baseline  for 
assessing  the  performance  of  the  near-neighbor  techniques.  An  interesting  realization  is 
that,  although  the  NN  algorithms  permit  actual  access  to  be  superior,  the  overall  cost  of 
the  NN  approach  could  be  greater  than  brute  force.  This  occurs,  for  example,  if  the  NN 
data  structure  must  be  set  up  too  frequently.  If  the  same  data  structure  can  be  used  for 
groups  of  contact  reports,  then  the  NN  algorithms  should  be  much  faster. 

Some  near-neighbor  algorithms  are  also  imperfect  in  that  some  near-neighbors,  by 
whatever  definition  of  “nearness,”  might  be  missed  if  the  search  is  too  restricted.  Clustering 
algorithms  can  permit  ambiguous  assignments  of  objects  to  particular  clusters,  for  example. 
The  MLG  algorithm  can  require  searches  of  irregular  regions  of  index  space  or  of  large 
regions  in  order  to  insure  finding  all  of  the  near  neighbors.  Tree-like  searches  of  the  MLG 
data  structure  [6]  might  be  necessary  in  some  cases. 

We  have  described  the  manner  in  which  NN  algorithms  can  be  used  in  correlating 
bearing  data  from  a  given  sensor  with  a  group  of  three-dimensional  tracks.  This  is  straight¬ 
forward,  and  appropriate  test  results  on  data  access  in  this  situation  are  available1  for  the 
MLG  and  Hockney’s  method. 

The  paper  has  also  presented  an  example  of  a  high-speed  algorithm  for  converting 
multisensor  bearing  data  to  positions  for  the  objects  within  the  region  common  to  all  of 
the  sensors.  These  algorithms  will  be  superior  to  any  brute  force  method  for  three  or 
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more  sensors,  since  the  latter  scales  as  Nm,  where  m  is  the  number  of  sensors,  while  our 
algorithm  is  particularly  suited  for  NN  algorithms,  which  scale  as  N  In  N.  The  algorithm 
also  reduces  the  problem  from  determining  the  m— fold  intersections  of  lines  to  the  pairwise 
equality  of  special  angles.  Thus  the  complexity  is  greatly  reduced.  Performance  tests  and 
estimates  indicate  that  the  algorithms  might  already  give  existing  hardware  the  capability 
to  satisfy  the  estimated  processing  requirements  for  situations  in  which  large  numbers  of 
objects  axe  moving  at  high  speeds.  Certainly  this  will  be  the  case  in  the  future,  given  the 
anticipated  improvements  in  speed  and  memory  of  computer  hardware. 
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Appendix:  Details  of  Multisensor 
Bearing-to-Position  Conversion 


For  completeness,  this  appendix  gives  details  on  the  transformation  of  the  bearing 
data  from  each  sensor  to  new  angles  denoted  by  $.  Figure  5  returns  to  the  example  used 
previously  (Fig.  3)  for  a  bearing  vector  ejy  corresponding  to  Sensor  B.  Figure  5b  shows 
that  the  angle  'c(j )  can  be  computed  from  the  unit  normals  to  the  plane  of  the  sensors 
and  the  plane  defined  by  et,c  (the  vector  from  Sensor  B  to  Sensor  C)  and  the  bearing  vector 
ebj- 

First  define  the  unit  vector  03  perpendicular  to  the  plane  of  the  sensors.  If  the  site 
producing  the  bearing  vector  is  l  S  (A,  B,  C},  then  we  define  the  values  of  l  —  1  and  l  4- 1 
by  the  preceding  and  following  sensor  indices  in  the  appropriate  cyclic  permutation  of  the 
sequence  A  —  B  —  C.  For  example,  if  /  =  B,  then  /  —  1  =  A  and  1  +  1  =  C.  Then  the 
desired  unit  vector  is  the  vector  cross  product 


0  =  e*«+i)  *  e(t~i)' 
l«i(f+i)  x  e(/-i)(l 

In  our  example  with  l  =  B,  this  becomes 


&bc  X  6a& 
|®bc  X  Ca&| 


The  perpendicular  to  the  plane  BCj  is 


Pbc(j) 


Bbj  x  e&c 
|®6/  X  e&c| 


The  angle  then  satisfies  the  equations 


(A.l) 


(-4.2) 


(A.3) 


sin  $6*c U)  =  (AcO')  x  /?,)  •  eic  (A. 4a) 

cos  =  Pbc (j)  ■  0,  ,  (.4.46) 
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in  which  the  centered  dot  (•)  signifies  a  scalar  product.  Figure  5c  shows  the  definition  of 
the  second  angle  corresponding  to  the  bearing  vector  e&y.  This  is  the  angle  between 

the  plane  of  the  sensors  and  the  plane  defined  by  the  line  joining  sensors  A  and  B  and  the 
bearing  vector  ejy.  By  convention  (Section  IV.  1),  the  observer’s  line  of  sight  runs  from 
sensor  A  toward  sensor  B,  producing  the  view  shown  in  Fig.  5c.  A  cyclic  permutation  of 
the  subscripts  A  —  B  —  C  will  define  the  angles  corresponding  to  bearing  vectors  produced 
by  sensors  A  and  B. 

The  remaining  question  is  that  of  transforming  the  uncertainty  (A#,  Ao)  to  uncer¬ 
tainties  in  the  new  angles  $.  Again  consider  the  angle  ^c(j)  the  bearing  vector 
ejy  =  (Qb(j),  4>b(j)).  Define  the  right-handed  coordinate  axes  (visualize  with  Fig. 5b) 

x  =  eic 

!/  =  e.j  ■  (.4.5) 

z&0. 

Then,  if  the  subscripts  on  the  bearing  angles  are  dropped  for  convenience,  the  important 
unit  vectors  axe 

0s 

e&c 

Equation(A.3)  then  gives  us 

0bc{j )  —  G  (0,  cos  9,  —  sin  9  sin  6)  ,  (-4.7) 

where 

G  =  (cos2  9  -f  sin2  <?sin2  <p)~ »  .  (-4.S) 


=  (sin0cos<£,  sin  9  sin  6,  cos9) 

=  (0,0,1)  .  (.4.6) 

=  (1,0,0) 


Equations(A.4)  then  become 
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sin^J^Q')  =  G  cos  9 
cos  =  —  G  sin  9  sin  <j> 


(.4.0) 


Equations  (A.7)-(A.9)  show  that  the  error  in  the  new  angles  varies  with  the  bearing  vector 
ejj.  The  new  angles  themselves  can  become  undefined  for  pathological  cases  in  which  the 
three  bearing  vectors  are  not  distinct  or  for  which  the  bearing  vectors  lie  close  co  the  plane 
of  the  sensors.  An  example  would  be  an  object  lying  along  or  near  the  line  ej,c.  in  which 
case  (<9,  4>)  =  (f ,  0).  At  the  other  extreme,  for  9  =  mir,  where  m  is  an  integer,  and  variable 
</>,  the  error  in  $^>0  )  is  A<52  ~  A02.  The  general  equations  for  the  error  become  unwieldy 
because  of  the  normalization  factor  G.  A  detailed  examination  is  necessary  to  determine 
the  ranges  of  bearing  vectors  for  which  their  intersections  yield  useful  position  data.  Since 
our  method  is  equivalent  to  the  brute  force  (“intersection'’)  approach,  the  range  of  validity 
will  be  the  same.  One  method  of  determining  this  range  in  the  context  of  our  method  is 
to  use  the  equations  in  this  appendix  in  conjunction  with  a  limit  on  the  allowable  size  of 
the  errors  in  the  angles  $.  An  example  of  such  a  constraint  is 


(A $scC/'))2  <  (A#  +  Ac?)2  . 


(.4.10) 


